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ABSTRACT 

We have analyzed new and archival time series spectra taken six years apart during transits of the 
hot Jupiter WASP-33 b, and spectroscopically resolved the line profile perturbation caused by the 
Rossiter-McLaughlin effect. The motion of this line profile perturbation is determined by the path 
of the planet across the stellar disk, which we show to have changed between the two epochs due to 
nodal precession of the planetary orbit. We measured rates of change of the impact parameter and 
the sky-projected spin-orbit misalignment of dh/dt = —0.0228^0 oois dX/dt = —0.487^0 076 ° 

yr“^, respectively, corresponding to a rate of nodal precession of dQ/dt = 0.373^0 083 ° This is 

only the second measurement of nodal precession for a confirmed exoplanet transiting a single star. 
Einally, we used the rate of precession to set limits on the stellar gravitational quadrupole moment of 
9.4 X 10-^ < J 2 < 6.1 X 10-^. 

Subject headings: line: profiles — planetary systems — planets and satellites: individual: WASP-33 b 
— planet-star interactions — techniques: spectroscopic 


1. INTRODUCTION 

WASP-33 b is a hot Jupiter orbiting a relatively mas¬ 
sive (1.5Mq), rapidly rotating {vsini^ = 85.6 km s“^) 
star (Collier Cameron et al. 2010b), which is notable 
for being one of the hottest planet host stars known 
(Teff = 7430 K). Due to the wide, rotationally broad¬ 
ened stellar lines. Collier Cameron et al. (2010b) were 
only able to set an upper limit on the radial velocity re¬ 
flex motion of the host star due to the planet. Even in 
the Kepler era, detection of this motion is typically nec¬ 
essary to confirm a transiting giant planet candidate as a 
bona fide planet. Instead, they confirmed the planetary 
nature of the transiting companion using Doppler tomog¬ 
raphy. This method relies upon the Rossiter-McLaughlin 
effect (Rossiter 1924; McLaughlin 1924), where the tran¬ 
sit of a companion across a rotating star causes a pertur¬ 
bation to the rotationally broadened stellar line profile. 
Doppler tomographic observations allow us to resolve this 
line profile perturbation spectroscopically, unlike typical 
radial velocity Rossiter-McLaughlin observations, where 
the perturbation is interpreted as an anomalous radial 
velocity shift due to the changing photocenter of the line 
(e.g., Triaud et al. 2010). Most importantly for this work, 
the movement of the line profile perturbation across the 
line profile during the transit maps directly to the path of 
the planet across the stellar disk, allowing us to measure 
the location and orientation of the transit chord relative 
to the projected stellar rotation axis. 

Using their Doppler tomographic observations of 
WASP-33 b. Collier Cameron et al. (2010b) measured a 
sky-projected spin-orbit misalignment of A = —105.8° ± 
1.2° (using their dataset from McDonald Observatory). 
They also found an orbital period of P = 1.2198669 ± 
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0.0000012 days. Since WASP-33 b is on a highly in¬ 
clined, short-period orbit about a rapidly rotating (and 
therefore likely dynamically oblate) star, lorio (2011) es¬ 
timated that the orbital nodes should precess at a rate 
of dVt/dt < 8.2 X 10“^^ s“^ (< 1.5° yr“^). They pre¬ 
dicted that this would result in a changing transit du¬ 
ration that would be detectable in ^ 10 years. Such 
a measurement will be challenging, however, as WASP- 
33 is a J Set variable (Herrero et al. 2011). The stel¬ 
lar non-radial pulsations cause distortions in the transit 
light curve, which could induce systematic errors in the 
measurement of the transit duration. This change in 
the transit duration, however, is caused by the changing 
impact parameter (denoted 6), which can be more accu¬ 
rately measured using Doppler tomography than using 
the transit lightcurve. Eor instance. Collier Cameron et 
al. (2010b) measured the impact parameter of WASP- 
33 b to be 6 = 0.176 ± 0.010 using their spectroscopic 
data and b = 0.155lo.i20 using their photometric data. 

It has now been more than six years since the Doppler 
tomographic observations presented by Collier Cameron 
et al. (2010b) were obtained. This offers a sufficient time 
baseline to allow the detection of the movement of the 
transit chord due to nodal precession. We have thus 
collected a second epoch of Doppler tomographic obser¬ 
vations, and have detected the changing transit chord. 

Orbital precession has previously been detected for 
Kepler-13 Ab by Szabo et al. (2012), who measured 
a rate of change of the impact parameter of db/dt = 
—0.016 ± 0.004 yr“^ using the changing transit duration 
in Kepler photometry. Like WASP-33 b, Kepler-13 Ab 
is a hot Jupiter orbiting a rapidly rotating star on an 
inclined orbit (Johnson et al. 2014). Barnes et al. (2013) 
proposed a large rate of nodal precession for the young 
hot Jupiter candidate PTEO 8-8695 b (van Eyken et al. 
2012), but this planet candidate is still unconfirmed (Cia- 
rdi et al. 2015). Our measurement of the orbital preces¬ 
sion of WASP-33 b is thus the second such measurement 
for a confirmed exoplanet orbiting a single star. 
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2. OBSERVATIONS AND METHODOLOGY 
2.1. Spectroscopic Observations and Analysis 

Collier Cameron et al. (2010b) observed one transit of 
WASP-33 b with the 2.7m Harlan J. Smith Telescope 
(HJST) at McDonald Observatory on 2008 November 12 
UT, and we reanalyzed these data. They also observed 
two other transits with other facilities, but we did not 
reanylize these data. We observed a second transit with 
the HJST on 2014 October 4 UT, 2,152 days (1,764 plan¬ 
etary orbits) after the first dataset. We obtained both 
datasets using the Robert G. Tull Coude Spectrograph 
(TS23; Tull et al. 1995), with a spectral resolving power 
of R = 60,000. The exposure length was 900 seconds for 
both datasets. We obtained 13 spectra in 2008 and 21 in 
2014; 10 spectra in each dataset were taken during the 
transit. 

Our methodology for preparing the time series spectra 
for Doppler tomographic analysis was substantially the 
same as that described in Johnson et al. (2014). We ex¬ 
tracted an average line profile from each spectrum using 
least squares deconvolution (Donati et al. 1997), and sub¬ 
tracted the average out-of-transit line profile to produce 
time series line profile residuals, which are most useful 
for analysis. Thanks to WASP-33’s brightness (U = 8.3) 
we were able to obtain very high quality line profiles; the 
standard deviation of the continuum was 0.010 of the 
depth of the line profile for the 2008 observations and 
0.0078 for the 2014 dataset. We modeled the line profile 
perturbation due to the transiting planet by numerically 
integrating the line profile from each surface element on 
the star over the visible stellar disk, accounting for the 
finite exposure time. 

The stellar non-radial pulsations caused a pattern of 
striations in the time series line profile residuals, compli¬ 
cating the analysis. In order to minimize this effect we 
exploited the fact that the pulsations propagate in the 
prograde direction, whereas the planetary orbit is ret¬ 
rograde (|A| > 90°). The frequency components due to 
the pulsations and the planetary transit thus tended to 
be separated in the two-dimensional Fourier transform of 
the time series line profile residuals. We constructed a 
Fourier filter by multiplying each complex element of the 
Fourier transform by unity if that element was in a region 
where there was power only from the transit signature, 
and zero if the element was in a region with significant 
power from the pulsations, with a Hann function tran¬ 
sition between the two regimes. We then performed an 
inverse Fourier transform on the filtered Fourier spec¬ 
trum. This successfully removed most of the effects of 
the pulsations. For best results we had to filter out low- 
frequency modes where there was power from both the 
pulsations and the transit; however, the high-frequency 
components were sufficient to reconstruct most of the 
transit signature. 

We obtained best-fit values of the transit parameters 
by exploring the likelihood space of model fits to the data 
using a Markov chain Monte Carlo (MCMC) with affine- 
invariant ensemble samplers (Goodman & Weare 2010), 
as implemented in emcee (Foreman-Mackey et al. 2013). 
We performed a joint fit to the time series line profile 
residuals from both 2008 and 2014, as well as a single 
spectral line, the latter in order to measure the vsini^ 
of WASP-33. For the single line we fit a rotationally 


broadened line profile to the Ba ll line at 6141.7 A, cho¬ 
sen because it is deep but unblended and unsaturated. In 
order to minimize the impact of the line profile variations 
we stacked all of our 2008 spectra. We did not utilize the 
2014 dataset because we could not obtain a good fit to 
the Ba ll line. For the time series line profile residu¬ 
als we fit a transit model computed as described above 
and passed through our Fourier filter. The MCMC had 
sixteen parameters: A, 5, and the transit epoch Tq at 
the two epochs, vsini^, Rp/R^, a/R^^ P, four quadratic 
limb darkening parameters (two each for the single-line 
and the time series line profile residual data), the width 
of the Gaussian line profile from an individual surface 
element (due to intrinsic broadening, thermal broaden¬ 
ing, and microturbulence), and a velocity offset between 
the single line data and the rest frame. All parameters 
except To, A and b were assumed to remain constant be¬ 
tween 2008 and 2014. We calculated Tq for 2008 from the 
epoch and period given by Kovacs et al. (2013), while the 
2014 transit epoch was taken from our simultaneous pho¬ 
tometric observations of that transit (see §2.2). Rather 
than fitting the limb darkening coefficients directly, we 
used the triangular sampling method of Kipping (2013). 
We set Gaussian priors upon P^/P^, a/P*, P, and the 
2008 To, and set the prior value and width to the parame¬ 
ter value and uncertainty, respectively, found by Kovacs 
et al. (2013), while the prior value and width for the 
2014 To were taken from our photometric observations. 
For the limb darkening parameters and prior values we 
used the same methodology as Johnson et al. (2014). We 
used a set of 100 walkers, ran each one for 1000 steps, 
and cut off the first 500 steps of convergence and burn-in, 
resulting in 50,000 samples from the posterior distribu¬ 
tion. 

2.2. Photometric Observations and Analysis 

Deviation of the observed transit midpoint from that 
expected based on the published ephemeris could mas¬ 
querade as a change of the transit chord in the Doppler 
tomographic data. With exposure lengths of 900 sec¬ 
onds, the spectroscopic data alone did not sufficiently 
constrain the transit epoch. In order to better constrain 
this parameter we simultaneously obtained photometry 
of WASP-33 using the Las Gumbres Observatory Global 
Telescope Network (LCOGT; Brown et al. 2013) Im tele¬ 
scope and SBIG camera at McDonald Observatory. 

We observed in the Sloan T band, and defocused the 
telescope in order to reduce the effects of inter-pixel vari¬ 
ations and avoid saturating on this bright star (V = 8.3). 
We obtained 700 images, each with an exposure length of 
10 seconds. We used the astrometry.net code (Lang et al. 
2010) to register all images, and then performed aperture 
photometry on WASP-33 and three reference stars using 
the IDL task APER. We calculated the formal uncer¬ 
tainty on each data point incorporating the uncertainty 
from APER (based upon photon-counting noise and un¬ 
certainty in the sky background), as well as an estimated 
contribution from the calibration frames and from scin¬ 
tillation noise (Young 1967). 

We produced a model of the transit lightcurve using 
the JKTEBOP package^ (e.g., Southworth 2011). We 
used an MCMC to produce posterior probability distri- 

^ http://www.astro.keele.ac.uk/jkt/codes/jktebop.html 
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butions for each of the model parameters (we used a cus¬ 
tom MCMC derived from that in Johnson et ah 2014, 
not that built into JKTEBOP). 

There were significant distortions of the lightcurve due 
to the aforementioned stellar variability. We treated the 
stellar variability as correlated noise, and modeled it non- 
parametrically using Gaussian process regression (e.g., 
Gibson et al. 2012 ; Roberts et al. 2013). This method¬ 
ology has been used previously to treat stellar noise in 
transit lightcurves (e.g., Barclay et al. 2015). 

We constructed the covariance matrix K using a 
Matern 3/2 kernel, where each element of the matrix 
was 


kij = 


1 + 


V^\u - tj 
i 


exp - 


V^\u — tj 
i 


+ ^ij^i 

( 1 ) 

where i^j denote two of the photometric observations, 
ti, tj are the times at which observations i,j were ob¬ 
tained, ai is the formal uncertainty on datapoint i, and 
a and I are hyperparameters describing the amplitude 
and timescale of the stellar variability, respectively. 

The MCMG used eight parameters: 6 , Rp/R^, ajR^^ 
(a, /, the epoch of transit center Tq, and two quadratic 
limb darkening parameters. We set Gaussian priors upon 
RpjR^ and ajR^^ using the best-fit values and uncertain¬ 
ties found by Kovacs et al. (2013) as the center and width 
of the priors, respectively. We also set priors upon the 
limb-darkening parameters, using JKTLD^ to find the 
expected limb darkening values from ATLAS model at¬ 
mospheres from Glaret (2004) at the stellar parameters 
of WASP-33 found by Gollier Gameron et al. ( 2010 b). 


3. RESULTS 

The LGOGT lightcurve is shown in Fig. 1. We found 
a best-fit time of transit center of Tq = 2456934.77146 ± 
0.00059 BJD. This is 12.3 minutes later than predicted 
by the ephemeris of Gollier Gameron et al. ( 2010 b), but 
is in agreement with that predicted by the ephemeris of 
Kovacs et al. (2013). 

We show the time series line profile residuals in Fig. 2 . 
The best-fit values of the model parameters are given 
in Table 1. We found that both the impact parame¬ 
ter and the spin-orbit misalignment have changed be¬ 
tween the two epochs: we measured b = 0 . 218 lo!o 29 
A = -llO.Oetn? ° in 2008, and b = 0.0840t^;l^^i^ and 
A = —112.93^0 2 ? ° in 2014. Our uncertainties on these 
values are rather small (cf. Gollier Gameron et al. 2010 b, 
whose uncertainties on A and b are ^2 — 5 times the size 
of ours). They did not remove the stellar pulsations from 
their data, and so it is perhaps not unexpected that we 
can obtain a more precise result. There may, however, 
be sources of systematic errors which were not taken into 
account in our calculation of the uncertainties. The val¬ 
ues of A and b that we obtained from the 2008 data dis¬ 
agree with those found by Gollier Gameron et al. ( 2010 b) 
by 3.4cr and 1.3 <j, respectively, but agree with another of 
their Doppler tomographic datasets to within 1.2cr. Most 
likely this results from differing treatments of the stellar 
pulsations. Additionally, our posterior distribution for b 
in 2008 is double-peaked, resulting in asymmetric uncer- 


^ http://www.astro.keele.ac.uk/jkt/codes/jktld.html 
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Figure 1. Top: LCOGT lightcurve of the transit of WASP-33 b 
on 2014 October 4 UT. The data are shown in gray, with the best-fit 
transit model in red and the best-fit transit plus Gaussian process 
model in blue. Middle: residuals with the best-fit transit model 
subtracted, showing the stellar variability. Bottom: residuals with 
the best-fit transit plus Gaussian process model subtracted, show¬ 
ing the power of the Gaussian process to model and remove the 
stellar variability. 


tainties on this parameter. 

Although we did not measure the transit duration di¬ 
rectly from our data, we calculated the expected duration 
using Eqn. 3 of Seager & Mallen-Ornelas (2003); this is 
shown in Table 1. The transit duration implied by our 
Doppler tomographic measurements has changed by 2.7 
minutes between the two epochs, a challenging measure¬ 
ment for typical ground-based data, even without the 
complication of stellar variability. 

Using our values of b and A at the two epochs, we calcu¬ 
lated the rate of precession. We used the definition of the 
argument of the ascending node U as given by Queloz et 
al. ( 2000 ), i.e., the angle between the plane of the sky and 
the intersection between the planetary orbital plane and 
a plane parallel to the line of sight which is also perpen¬ 
dicular to the projection of the stellar rotation axis onto 
the plane of the sky, as measured in this latter plane. See 
Fig. 2 of Queloz et al. (2000) for a graphical definition; 
note that the quantity they denoted as A is our 6 , and 
their i is our A- Using this definition and the definition 
of the impact parameter, b = a/R^cosip, we related U 
to our known quantities with 

tan U = — sin A tan ip ( 2 ) 

which we used to calculate U at the two epochs. We 
assumed that a/R^ remains constant. We found db/dt = 
-0.0228to;“i8 yr“i and dX/dt = -0.487to;o76 ° yr“\ 
and calculated a rate of nodal precession of WASP-33 b 
of dCl/dt = 0.373lo.o83 ° This is in agreement with 

the prediction of lorio (2011), dCl/dt < 1.5° yr“^. A 
schematic view of the changing transit chord is shown in 
Fig. 3. 

The observation that both b and A are changing implies 
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Figure 2. Doppler tomographic datasets and Fourier filters. The left and right columns show the 2008 and 2014 datasets, respectively. 
Top row: raw observations. Second row: the data after application of the Fourier filter. Each colorscale line denotes the deviation of the 
line profile at that time from the average out of transit line profile; bright areas denote shallower areas of the line. Time increases from 
bottom to top, and we define the “transit phase” such that it equals 0 at ingress and 1 at egress. Vertical dashed lines mark r? = 0, sinz*, 
a horizontal dashed line marks the time of mid-transit, and the four small crosses denote the times of first through fourth contacts. The 
transit signature is the bright streak running from bottom center towards the upper left. The time range depicted is the same for all 
plots; flat blue areas indicate regions where we do not have any observations. Most of the remaining anomalous structure in the filtered 
datasets is ringing due to the filter. The transit signature has shifted slightly to the right between the 2008 and 2014 epochs. Third 
row: two-dimensional Fourier transforms of the time series line profile residuals, shown with a square-root color scale to best display the 
frequency structure. The transit signature is the narrow structure running from upper left to lower right. Bottom row: masks used to 
Fourier filter the data. 


that the total angular momentum vector of the system 
Ltot (the sum of the stellar spin and planetary orbital 
momentum vectors), about which the planetary orbital 
angular momentum is precessing, is neither close to per¬ 
pendicular nor close to parallel to the line of sight. Con¬ 
sider two limiting cases: if Ltot were perpendicular to the 
line of sight, then the precession would manifest as purely 
a change in 6 , whereas if Ltot were parallel to the line of 
sight, the precession would manifest as purely a change 
in A. Intermediate motion implies an intermediate angle. 

Using the rate of precession we set limits on the stellar 
gravitational quadrupole moment J 2 . This is 


__dn^ 

^ dt Sir 



(3) 


(e.g, from rearranging Eqn. 10 of Barnes et al. 2013), 
where ip is the angle between the stellar spin and plane¬ 
tary orbital angular momentum vectors (A is projected 
onto the plane of the sky). This angle can be expressed 
as (from Eqn. 25 of lorio 2011) 

cos pj = cos cos ip -h sin sin ip cos A (4) 

and, while we do not know U or following lorio ( 2011 ) 
we can set limits on these quantities. By requiring that 


WASP-33 rotate at less than the breakup velocity, and 
using the stellar parameters found by Collier Cameron 
et al. ( 2010 b), lorio ( 2011 ) set limits of 11 . 22 ° < < 

168.77°. Along with our values of ip and A, this implies 
a Icr range of 93.06° < < 110.33°. Thus, we set limits 

of 9.4 X 10“^ < J 2 < 6.1 X 10“^. Eor comparison, the 
Solar value is J 2 ^ 2 x 10“^ (e.g., Roxburgh 2001). 

4. CONCLUSIONS 

We have detected the nodal precession of the hot 
Jupiter WASP-33 b, and measured a rate of change of 
the ascending node of dU/dt = 0.373lo.o83 ° This 

implies that WASP-33 b began transiting its host star 
as viewed from the Earth in 1974^3, and will transit 
until 2062 ^^ 0 ; the precession period is ^ 970 years. Eur- 
thermore, we have set limits on the stellar gravitational 
quadrupole moment of 9.4 x 10“^ < J 2 < 6.1 x 10“^. 

Given the rate of change of the impact parameter of 
Kepler-13 Ab found by Szabo et al. ( 2012 ) and the uncer¬ 
tainty on b measured by Johnson et al. (2014), we expect 
that the precession of Kepler-13 Ab will be measurable 
by Doppler tomography by 2017. This will be important 
as Masuda (2015) found that such a measurement will be 
able to unambiguously distinguish between the conflict- 
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Table 1 

Observed Parameters 


Parameter 

2008 

2014 

V sin G (km s“^) 
To (BJD) 

86.63l°J^ 

2456934.77146 ± 0.00059 

a 


0.00173 ± 0.00082 

1 (minutes) 


20.7 ±9.2 

A(°) 


-112.93lo 21 

b 

0 91 Q+0-011 

u.^j.o_o 029 

0.0840lO;00?o 

ip (°) 

86.6llo 17 

88.695l°;“i 

D(°) 

86.39l“l8 

88.584l°;°34 

ri 4 (days) 

0 . 11694 lHooI? 

0.11880 ±0.00033 


Note. — Uncertainties are purely statistical and do not take 
into account systematic sources of error. The transit duration ri 4 is 
calculated from b, P, and a/R^, and is not measured directly. 



Figure 3. Schematic showing the transit chord crossing the star 
at the 2008 and 2014 epochs. The stellar rotation axis is vertical, 
and the north pole is at the top, such that star rotates from left to 
right. The planet moves along the red lines from bottom to top. 
The silhouette of the planet is shown for a 2008 transit. The limb 
darkening shown on the star corresponds to that for our best-fit 
Doppler tomographic solution. Surface brightness variations due 
to the stellar pulsations are not depicted. 

ing values of the system parameters found by Johnson 
et al. (2014) via Doppler tomography and Barnes et ah 
(2011) using the effects of stellar gravity darkening on the 
lightcurve. No other known planet is currently amenable 
to the detection of orbital precession with Doppler to¬ 
mography, as all other planets with published Doppler 
tomographic observations have approximately aligned or¬ 
bits and thus should display much slower precession than 
WASP-33 b (Collier Cameron et al. 2010a; Miller et al. 
2010; Gandolfi et al. 2012; Brown et al. 2012; Bieryla et 
al. 2015; Bourrier et al. 2015). Current and future tran¬ 
sit surveys can provide more targets amenable for the 
detection of precession via Doppler tomography over the 
next decade. 
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